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Abstract 

In this article, a formulation of a point-collocation method in which 
the unknown function is approximated using global expansion in tensor 
product Bernstein polynomial basis is presented. Bernstein polynomials 
used in this study are defined over general interval [a,b]. Method incor- 
porates several ideas that enable higher numerical efficiency compared to 
Bernstein polynomial methods that have been previously presented. The 
approach is illustrated by a solution of Poisson, Helmholtz and Bihar- 
monic equations with Dirichlet and Neumann type boundary conditions. 
Comparisons with analytical solutions are given to demonstrate the ac- 
curacy and convergence properties of the current procedure. The method 
is implemented in an open-source code, and a library for manipulation of 
Bernstein polynomials bernstein-poly, developed by the authors. 



1 Introduction 

Bernstein polynomials are known because of their many useful properties pQ. 
When restricted to the unit interval they are used in the definition of Bczicr 
curves, which are very important tools in computer graphics. On their own, 
they have enjoyed increased attention in recent years, specifically as means to 
represent solutions to differential equations. Bhatti and Bracken [5] present a 
Galerkin method which uses Bernstein polynomials as trial functions for solution 
of two-point boundary value problem, Yousefi et. al [3] use Bernstein polyno- 
mials in Ritz-Galerkin method to approximate solution of hyperbolic PDE with 
an integral condition. Doha et. al. 0] present the solution method for high- 
even-order ordinary differential equations. 

Unlike all listed examples, which use Bernstein polynomials in some sort of 
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Galerkin method, in this article we develope a collocation method. The mo- 
tivation to do so comes from the nature of Bernstein polynomials themselves. 
When function is expanded in Bernstein polynomial basis, the expansion coeffi- 
cients also respresent the nodal values of the expanded function. For two-point 
boundary value problems this leads to very simple implementation, allowing 
direct imposition of the boundary conditions at the end-points of the inter- 
val. When the unknown function is approximated as a global expansion in 
Bernstein polynomials basis, and used in a point-collocation method, the ap- 
proach becomes very similar to pseudospectral methods [5] , [6] [7] . Unlike their 
polynomial counterparts used in pseudospectral methods (e.g. Chebyshev and 
Legandre), Bernstein polynomials are not orthogonal, and many examples of 
basis transformation exist when this is necessary (eg. [S]). 
For the present article we assume that L, the elliptic differential operator in 
interior, and B, the boundary operator are linear, and that they define well 
posed boundary value problem 

L[u](x) = /(x), xe!l. (1) 

B[u](x) = .g(x), xecKl (2) 
Where f2 £ M 2 is a given rectangular domain. 

In the remainder of the paper we give a brief review of Bernstein polynomial 
properties, relevant to the proposed method. The solutions to PDEs defined 
over two-dimensional domains can be represented by surfaces embedded in 
three-dimensional Eucledian space. This serves as rationale of Section 3, in 
which we describe surfaces defined by expansions in tensor product Bernstein 
polynomial basis, also giving expressions for elliptic operators that we will use 
subsequently. Then, in Section 4, we describe the fomulation of the proposed 
collocation method. Finally we give couple of example solutions, as well as their 
error, and convergence analysis in Section 5. 



2 Properties of Bernstein Polynomials 

In further discussion we consider only generalized Bernstein polynomials, those 
defined over arbitrary interval [a, b], and simply call them Bernstein polynomials. 
The Bernstein polynomials of nth degree form a complete basis over [a, b], and 
are defined by 

/ \ fn\ (x — a) l (b — x) n — i , . 

B i , n (x) = [ i y , i = 0,l,...,n, (3) 

where binomial coefficients are given by 

fn\ _ nl 

i\(n-i)V 1 ' 

They satisfy symmetry Bi }U (x) = B n -i tTl (l — x), positivity Bi }TL (x) > 0, and 
they form partition of unity X)"=o Bi_ n {x) — 1 on defining interval x £ [a, b]. For 
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i = 0, and i = n, they have value equalt to one at x = a, and x = b, respectively. 
Otherwise, they have a unique local maximum occuring at x = i/n, having the 
value B it n(i/n) = i l n- n {n - 

In a recent article [3] , authors present explicit formula for calculation of arbitrary 
order derivatives of Bernstein basis functions of any order defined over standard 
interval [0,1]. In an atempt to formulate collocation method for solution of 
BVP's over arbitary intervals, first step would be to write down this expression 
for Bernstein polynomials defined for arbitrary interval. 

Let Bi tn (x) define i — th basis function of n — th order, at point x, where 
x G [a, b], and i € [0,n]. Then, the derivative of order p of such basis function 
can be expressed as: 

j min(i,p) . . 

DPB ^= { n-m-ay E 

v /v ' k=max(0A+p-n) v ' 

Useful proprerty of Bernstein basis functions is that they all vanish at end- 
points of the interval, except the first and the last one, which are equal to one 
at x = a and x — b respectively. 

The factorial formula Eq. [4] for binomial coefficients is known not to be most 
numerically efficient for large numbers, and the alternative representations exist, 
most numerically efficient being the multiplicative formula, which we will use 
here 



3 Surfaces defined by tensor product of Bern- 
stein polynomials 

Following equations will help us in developing a collocation method using Bern- 
stein polynomials for boundary value problems defined over two-dimensional 
domains. 

Let the following equation represent a surface in K 3 defined by a tensor 
product of Bernstein basis functions 

n m 

f(x,y) =J2J2^Bi,n^)B j ,m(y)- (6) 

i=0 j=Q 

We can define p-th order partial derivative in x -axis direction in a following 
way 

& ^%r 1 = E E falPBiMB^y), (7) 
j=o j=a 

and an q-th order derivative in y -axis direction 
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d q f(x,y) 



^^1,., //,„;,-)//' «,,„(,/;. (8) 



dyi 

y *=o j=0 

Finally mixed derivative of order p + q is defined as follows: 

dP+q d f J x J ] = E E o) 

y »=0 j=0 

This paper deals with elliptic equations, therefore we will define two differential 
operators, most common in second and forth order boundary value problems. 
Laplacian operator can be written as 

Af(r v) - + gVfag) _ 

/l,2/J " dx 2 + dy? ~ 

n m n rn 

= EEftA m w fl25 .,»( i )+EEftj^»w i)2 v. (io) 

which can further be simplified to 



n m 

Af(x,y) =EE 5 ^ [%m(2/)£> 2 -Bi,n(z) + B i!n (x)D 2 B j>m (y)] . (11) 

Bicharmonic operator or bilaplacian can be written as 

A 2,/ s _ d A fi x ,y) , ^f{x,y) d i f(x,y) _ 
I[ ,V) ~ dx* + dxidyi + 9y 4 ~ 

n m 

= E E fa [ B jMv)D 4 B i>n (x) + 2D 2 B l , n (x)D 2 B J , m (y) + B^ n {x)D A B^ m {y)) . 

j=0 3=0 

(12) 

4 Collocation method formulation 

In present method, it is assumed that a variable can be expressed as an approx- 
imation in the form of global expansion in tensor product Bernstein polynomial 
basis. 

We look for the approximate solution in the form |6| for (x, y) € ft = ft U dft. 
The expansion coefficients ftj are unknown, and need to be determined. 
The collocation method, that we use to find the unknown coefficients, is a nu- 
merical procedure in which we require that solution satisfies differential equation 
exactly in a discrete set of points, known as collocation points. Number of col- 
location points has to match the number of unknowns. 
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If the order of Bernstein polynomial is n and to in x and y-axis directions re- 
spectively, then there is n c = (n + 1) x (m + 1) unknown coefficients in global 
polynomial expansion. If we are solving homogenous Dirichlet problem the num- 
ber of unknown is n c d — (n — 1) x (to — 1) . Reason for this is a property of 
Bernstein polynomials that we take the advantage of, namely the property to 
vanish at end points, except for the first and the last one, that we mentioned in 
Section [2] If the function value at the boundary is equal to zero, we eliminate 
the non-vanishing basis functions by setting their associated expansion coeffi- 
cients to zero, and continue by solving the problem for the smaller number of 
unknowns n cc i. 

Discretized equation set is obtained by substituting the unknown function and 
it's derivatives with the appropriate representation in Bernstein polyniomial 
basis, as described in Section . Substituting these Bernstein polynomial inter- 
polants into the original equation set, and applying it at the set of collocation 
points, one arrives at a linear system of equations 

Ab = c. (13) 

An element of the system matrix a^j is obtained by evaluating the differen- 
tial operator on the tensor product of ith and jth Bernstein polynomial basis 
function, at the collocation point defined by Cartesian coordinate pair (xi,yt). 

a i,j = L [ B i,n{%i) B j,m{Vi)], i = 1, n - 1, j = 1, .., m - 1. (14) 

An element of the right-hand side vector Cj is an evaluation of the right- 
hand side function at the ith collocation point defined by Cartesian coordinates 

(xi,Vi). 

Ci = f(xi,Vi), i = l,...,n cd . (15) 

For the cases other than those with homogenous Dirichlet boundary condi- 
tions, the number of unknowns, as being said, is larger and linear system has 
the additional rows defined by 

dij = B [Bi, n (xi)Bj, m {yi)] , (16) 

Ci = g(xi,yi), (17) 

Where coordinates (xi,y{) define colocation points that belong to the domain 
boundary. Index pairs belong to the set (0, l...m — I) U (n, l...m — I) U 

(l...n - 1, 0) U (I. ..ri - 1, to) U (0, 0) U (n, 0) U (0, to) U (n, to). There is in total 
2(n + to) additional rows required to specify boundary conditions. 
The change in boundary conditions ([2|, amounts to changing a few rows in the 
matrix A, as well as in the right-hand side vector c. Using expressions defined 
in Section 3. we may define any type of boundary conditions. In particular, the 
operator B takes the form of Eq. (|6]) when we need to impose Dirichlet boundary 
conditions, or ([7|,([8]), with p = q = 1 for Neumann boundary conditions. 
Solution vector b is as an auxilliary one-parameter array, it's values are moved to 
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two-parameter array of Bernstein polynomial expansion coefficient /3. Mapping 
has the following form 

binode: (-^) 

Where inode = (m + 1) * (i — 1) + j, and i = 0, = 0, ...,m. If only 

inner expansion coefficients are sought for, in the case of homogenous Dirichlet 
problems, coefficient mapping takes the same form (18), but the indexes are 
defined by inode = (m — 1) * (i — 1) + j — 1, and i — 1, n — 1, j = 1, m— 1. 

Finally we note that the distribution of points on a tensor product grid may 
be uniform, or non-uniform depending on situation. The full linear system [13] 
is solved using LU decomposition with partial pivoting. 



5 Example Application 

Collocation method formulated in previous section is implemented in bernstein- 
poly, a library for manipulation of Bernstein polynomials. It is an open-source 
code developed by the authors, written in Python [5J. In what follows, we will 
present couple of examples, with the purpose to illustrate validity and accuracy 
of the present collocation method. In all examples we use the same order of 
Bernstein polynomials in both x and y axis direction. 

Example 1. 

Consider Poisson equation 



Au(x,y)=f(x,y), x, y G [-1, 1], (19) 

with homogenous Dirichlet boundary conditions. Source term is defined as 
f(x,y) = — 2ir 2 sin(irx) sin(iry). Exact solution for the given problem is 

Uexact = sin(nx) sin{iry) (20) 

Fig. 1 shows numerical solution for approximation using Bernstein polynomials 
of order n=21. The absolute numerical solution error abserr(i, j) = \\u(i,j) — 
ufjH , where u\ ^ represents the exact solution, is estimated at each (i,j)th 
collocation node, and the result is presented in Fig. 2. To study how accuracy 
changes with incresing the order of polynomial approximation we use L 2 relative 
error norm, defined in the following way 

^rror= L '^; e ; j) . (21) 



Results of the analysis of L 2 relative error norm are sumarized in Table 1. 
To study order of accuracy it is useful to plot L 2 relative error norm as a 
function of polynomial order, n, which is done in log-log plot in Fig. [3] We 
observe exponential decay in the error for polynomial orders up to n — 17, after 
which the error decrease slows down. Eventually for order of Bernstein basis 
polynomials higher than n = 21, after which the error continually grows. 
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Figure 1: Bernstein polynomial collocation method solution of Poisson equation 
(Example 1), n=20. 



Example 2. 

In this example we consider the Poisson equation with non-homogenous Dirichlet 
boundary conditions. Treatment of this problem differs from the previous one 
algorithmically in a way we impose boundary conditions. For every collocation 
point at the boundary we write an additional linear equation and we solve linear 
system of size (n c ) x (n c ). 



Source term in the Poisson equation ( 19 1 for this problem is defined by f{x, y) = 
6xy(l — y) — 2a; 3 , and solution domain by unit square x,y G [0, 1]. The exact 
solution for this problem is 

Uexact = 2/(1 - y)x 3 (22) 

The figures |4j and [5] show numerical solution and absolute error distribution for 
the case of n = 20. very high accuracy is achieved quite "early", for n = 12 it 
reaches order of 10 , as seen in Table |J 

Example 3. 

Next example problem is defined by the Helmholtz equation 

(A + X)u = f(x,y). (23) 

The problem is originally found in [TU], and is defined by A = 1, f(x,y) = x, 
and non-homogenous Dirichlet boundary conditions which are derived from the 
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n 


L 2 rel. error 


11 


1.171 x 10 


-5 


13 


3.170 x 10" 


-7 


15 


6.536 x 10 


-9 


17 


1.049 x 10" 


10 


19 


5.137 x 10" 


11 


21 


3.111 x 10" 


11 


23 


3.966 x 10 


-9 


30 


3.146 x 10" 


08 


41 


3.499 x 10- 


07 


51 


1.091 x 10- 


05 


61 


1.140 x 10- 


04 


71 


1.346 x 10- 


04 



Table 1: L 2 relative error norm for the Example 1. 



n 


L 2 rel. error 


12 


5.841 x 10~ 15 


14 


2.595 x 10~ 14 


16 


1.754 x 10~ 13 


18 


5.039 x 10~ 12 


20 


1.544 x 10~ 10 


30 


1.907 x 10~ 8 



Table 2: L 2 relative error norm for the Example 2. 
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Figure 2: Absolute error (Example 1), n=20. 

exact solution. Domain of solution is square x, y € [— tt, ir\. Using Bernstein 
polynomials denned over general interval proves practical for this problem, as 
mapping of the domain to the unit square is not necessary. 
This problem is exactly solvable and the solution is 

Uexact = sin(x) + sin(y) + x. (24) 

As in previous examples we show numerical solution Fig. [6] and the distribution 
of the absolute error Fig. [7j Table [3] lists 1? relative error norm variation with 
increasing order of polynomial basis. 

Example 4. 

In next example we give the solution to biharmonic equation, which is often 
encountered in the theory of ellasticity, as it describes deflections of loaded 
plates. 

A 2 u(x,y) = f(x,y). (25) 

Two types of boundary conditions can be defined for this problem (I) u and 
d 2 u/dn 2 or (II) u and du/dn. Both cases are interesting on their own because 
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Figure 3: Example 1. L 2 relative error norm. 



n 


L 2 rel. error 


12 


9.035 x 10" 6 


14 


2.430 x 10~ 7 


16 


4.992 x 10~ 9 


18 


8.107 x 10~ n 


20 


4.057 x 10~ n 


22 


1.051 x 10~ 9 


30 


1.533 x 10~ 7 



Table 3: L 2 relative error norm for the Example 3. 




they require different solution approach. We will split discussion on biharmonic 
equation in two parts, present Example we will treat Type I and in the next 
one Type II boundary conditions. 

In the case of Type I boundary conditions, biharmonic equation can be split 
into two coupled Poisson equations 

Av(x,y) = f(x,y), (26) 

Au(x,y)=v(x,y). (27) 

The first example solution of biharmonic equation, taken from [TT], will deal with 
the case of simply supported rectangular plate < x < a,0 < x < b. Homoge- 
nous boundary conditions for both the function value and it's second derivatives 
in the direction normal to boundary are prescribed. 

Source function, that describes load distribution over the surface in theory of 
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0.8 



Figure 4: Solution of Poisson equation, non-homogenous Dirichlet BCs (Exam- 
ple 2), n=20. 

plates is given by 

a ( m 2 n 2 \ 2 mux nny 
f{x,y)=n I ^ + ^2 I sin — sin — . (28) 

Biharmonic equation defined in such a way, allows the exact solution 

mirx niry . . 

u exact (x,y) = sm—sm—r-. (29) 

For our purpose we set values m = n = l, a = b = l. Example solution for 
the order of polynomial, n = 20 is shown in Fig [8} The peaks in the absolute 
error near the corner nodes, as noticed in previous examples, is inherent to the 
present method, therefore we set an upper limit on the vertical axis here, to 
be able to get better picture of the distribution of absolute error in the rest of 
domain, Fig. [9] Table [4] shows variation in the error norm, in which the trend 
conforms to the one in previous examples. 

Example 5 

Biharmonic equation that is solved in this example has Type II boundary condi- 
tions and cannot be reformulated as a system of two coupled Poisson equations. 
This example is also taken from [TTJ. We will consider two cases, one having 
the exact solution 
Case 5a: 

f(x, y) = 56400(a 2 - Wax + I5x 2 )(b - y) 2 y A 
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1.05738e-10 - 
8.45902e-ll 
6.34427e-ll 
4,22951e-ll 
2.11476e-ll 




0.2 



0.4 



Figure 5: Absolute error (Example 2), n=20. 



n 


L 2 rel. error 


10 


6.538 x 10" s 


12 


4.854 x 10~ 10 


14 


2.832 x 10~ 12 


16 


1.284 x 10~ 12 


20 


9.467 x 10" 11 


30 


5.939 x 10~ 8 



Table 4: L 2 relative error norm for the Example 4. 



+ 18800a; 2 (6a 2 - 20a.x + 15.x 2 )y 2 (6fe 2 - 206y + 15y 2 ) 

+ 56400(a - x) 2 x i {b 2 - 10&y + 15y 2 ), (30) 
And one where the exact solution is unknown Case 5b: 

f(x, y) = Po = const, (31) 

where po takes value of 1000. The notation has physical significance, Example 
5b represents the case of a plate clamped at all four sides, and exposed to 
uniform load of fluid pressure. Example 5a admits the exact solution 

u eX act(x,y) = 2350x i {x~ a) 2 y\y-b) 2 . (32) 

In the assembling procedure, instead of writing the Eq. [l4]for the nodes laying 
right next to the egde nodes, we use Eq. |16[ which to remind once again, uses 
(x, y) values of the points on the boundary edge. These edge points are all 
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Figure 6: Solution of Hclmholtz equation, non-homogenous Dirichlet BCs (Ex- 
ample 3), n=20. 




- 3 14 -3.14 



Figure 7: Absolute error (Example 3), n=20. 
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Figure 8: Solution of Biharmonic equation (Example 4), Type I BCs, n=20. 




Figure 9: Absolute error (Example 4), n=20. 



14 



0.8 



Figure 10: Solution of Biharmonic equation (Example 5a), Type II BCs, n=20. 



11 


L 2 rel. error 


8 


2.506 x 10" 14 


10 


1.064 x 10" 14 


12 


3.773 x 10- 13 


14 


1.174 x 10~ n 


20 


3.853 x 10- 10 



Table 5: L 2 relative error norm for the Example 5a. 



immediate neighbours of those collocations points, we would normally write the 
equation for. We note that the Neumann boundary conditions are not written 
for the corner nodes. Numerical results, consistent with previous examples, is 
shown in Fig. [TUJ Fig. [IT] and Table [5] When there is no exact solution, we need 
to set up a criteria what solution to except as converged one. In all previous 
examples, which allowed exact solution, we see that the high-order accuracy is 
achieved by a comparatively small number of nodal points. For 21 nodal point 
in each direction (order of polynomial n=20), the L 2 relative error norm is of the 
order 10 10 to 10 11 . We should have the additional confidence in the result, if the 
spatial variation of solution expressed in terms of local maxima and minima, 



is not significant within the domain. Fig. 12 shows numerical solution with 



Bernstein polynomials or order n=20 in each direction. 
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Figure 11: Absolute error (Example 5a), n=20. 




Figure 12: Solution of Biliarmonic equation (Example 5b), Type II BCs, n=20. 
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6 Conclusions 



In this paper a novel formulation of the collocation method using Bernstein 
polynomials is proposed. The main reason why the collocation method is cho- 
sen, are it's flexibility and simple implementation. The methodology presented 
in this paper has been implemented in bernstein-poly code [SJ, and several ex- 
amples have been shown, where the elliptic boundary value problems in two 
dimensional domains were succesfully solved. Numerical results obtained by 
this method were compared with existing exact solutions. Excellent agreement 
and high accuracy is achievent even with small number of basis polynomials. 
Trough extstensive testing we have concluded that the three components of 
the algorithm used here: defining polynomials over general interval, the non- 
recursive formulation for derivatives and the use of multiplicative formula for 
individual binomial coefficients significantly enhance capabilities of the present 
procedure related to the previous Bernstein polynomial methods in terms of the 
flexibility and speed. 

As we have seen in log-log plots of L 2 relative error norm in all examples, the 
error decreases exponentially as the order of polynomial increases, until, around 
n = 17 — 20, when it changes the character and continually increases for higher 
values of n. We noticed the same character of error variation with n in all 
the examples above, which is related to the loss in accuracy in floating-point 
arithmetics with large numbers originating from factorials in the definition of 
Bernstein polynomials. This suggest further direction of investigation - devel- 
oping a method based on a principle of domain decomposition. In that case in 
each "element" we would keep moderate degree in Bernstein polynomial, and 
the method would resemble spectral element method. Similar has been done 
e.q. [12] . where Chebyshev multidomain method has been put forward. That 
study has served as a basis for development Spectral Difference Method. An- 
other advantage of the domain decomosition approach is that it produces sparse 
matrices, with the sparsity pattern usual for tensor product grid discretizations. 
Future plan is to develop bernstein-poly to the point where highly accurate so- 
lution of Navier-Stokes equations in complex three-dimensional domains is pos- 
sible. 
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Appendix A. Implementation of a new BVP in 
bernstein-poly 

Setting up a new elliptic boundary value problem in bernstein-poly code [9] is 
briefly described here on the case of Poisson equation presented in the Example 
3. Main program is located in collocation_test_2D.py, where the user needs to 
define solution domain by specifying it's x and y axis extents, and the order of 
approximating polynomials - n. 

Listing 1: Setting up domain borders and approximating polynomial order. 

# Bernstein polynomial order — the same order in both directions. 
n = 12 

m = 12 

# The number of unknowns — all nodes counted for the case with non— homogenous BCs. 
nvar = (n+l)*(m+l) 

# Solution interval [xl,x2] x [yl,y2]. 
xl = 0. 

x2 = 1. 

yl = 0. 
y2 = 1. 

# Uniform mesh. 

nd = Iinspace(xl,x2,n+1) 

Specifying the problem generaly described by Eq. (JlJ is straight-forward: 
Listing 2: Defining RHS vector entry 

def rhs(x.y): 

return — (6*x*y*(l— y)— 2*x**3) 



Listing 3: Defining system matrix entry 

def Ihs(i,j,n,m,xl,x2,yl,y2,x 1 y): 

return — Iaplacian(i 1 j,n,m,xl,x2,yl,y2,x 1 y) 



Listing 4: Defining the exact solution. 

def exact_solution(x,y): 

return y*(l— y)*x**3 

Boundary conditions have to be specified next. First we'll have a look how 
does the LHS matrix and RHS vector assembly look like (Listings 5 and 6). 

Listing 5: Forming RHS vector in the main function. 

f = zeros(nvar) # RHS vector 
for i in range(0,n+l): 
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3 x = nd[i] 

4 for j in range(0,m+l): 

5 y = ndfj] 

6 node = (m+l)*(i— l)+j 

7 # Non— homogenous BCs... 
s if(i==0): 

9 # Left side:: Run trough all betas 

10 f[node] = bc_left(x,y) 

11 elif (j==0): 

12 # Bottom:: Run trough all betas 

13 f[node] = bc_bottom(x,y) 

14 elif (i==n): 

15 # Right: Run trough all betas 

16 f[node] = bc_right(x,y) 

17 elif (j==n): 

is # Top: Run trough all betas 

19 f[node] = bc_top(x,y) 

20 else: 

21 f[node] = rhs(x.y) 



Listing 6: Forming system matrix 

1 K = zeros( (nvar.nvar) ) # LHS matrix 

2 for i in range(0,n+l): 

3 x = nd[i] 

4 for j in range(0,m+l): 
s y = ndfj] 

e node = (m + l)*(i— l)+j # node defines specific location on a grid 

7 # Non— homogenous BCs... 

8 if(i==0): 

9 # Left side:: Run trough all betas 

10 for k in range(0,n+l): 

11 for I in range(0,m + l): 

12 jfun = (m+l)*(k-l)+l 

13 Kfnodejfun] = Dirichlet(l,k,n,m,xl,x2,yl,y2,x,y) 

14 elif (j==0): 

is # Bottom:: Run trough all betas 

16 for k in range(0,n + l): 

17 for I in range(0,m + l): 

is jfun = (m+l)*(k-l)+l 

19 Kfnodejfun] = Dirichlet(l,k,n,m,xl,x2,yl,y2,x,y) 

20 elif (i==n): 

21 # Right: Run trough all betas 

22 for k in range(0,n + l): 

23 for I in range(0,m + l): 

24 jfun = (m+l)*(k-l)+l 

25 Kfnodejfun] = Dirichlet(l,k,n,m,xl,x2,yl,y2,x I y) 

26 elif (j==n): 

27 # Top: Run trough all betas 

28 for k in range(0,n + l): 



19 



for I in range(0,m+l): 

jfun = (m+l)*(k-l)+l 

K[node,jfun] = Dirichlet(l,k,n,m,xl,x2,yl I y2,x I y) 

else: 

# Interior: Run trough all betas 

for k in range(0,n+l): 

for I in range(0,m+l): 

jfun = (m+l)*(k-l)+l 

K[node,jfun] = Ihs(l,k,n,m,xl,x2 1 yl,y2 1 x,y) 



Listing 7: System matrix entry originating from the boundary operator ([2]) 

def Dirichlet(i,j I n,m,xl,x2 I yl,y2 I x,y): 

# Matrix entry for Dirichlet BCs: 

return basis_fun_eval(i,n,xl,x2,x)*basis_fun_eval(j,m,yl,y2,y) 

def Neumann_y_dir(i,j,n,m,xl,x2,yl,y2,x,y): 

# Matrix entry for Neumann BCs: 

return basis_fun_eval(i,n,xl,x2,x)*basis_fun_der(l I j,m I yl,y2,y) 

RHS vector entry for the collocation point located at domain boundary is 
evaluated according to function g{x) Eq. |2|. 

Listing 8: Evaluating RHS vector entry for the points on domain boundary, 
def bc_right(x,y): 

# Boundary condition for right edge of the rectangle 

return y*(l— y) 

The beauty of the proposed method and it's Python implementation, is that 
non-trivial problems descibed using PDEs are solved in a very small number of 
command lines. The guiding principle behind the development of the code is 
modularity and escalation towards solution of problems with greater complexity. 
Readers are encouraged to use and upgrade bernstein-poly in their own research. 

References 

[1] G. G. Lorentz. Bernstein Polynomials. University of Toronto Press, Toronto, 
Canada, 1953. 

[2] M. I. Bhatti and P. Bracken. Solutions of differential equations in a Bern- 
stein polynomial basis. Journal of Computational and Applied mathematics, 
205(l):272-280, 2007. 

[3] S.A. Yousefi, Z. Barikbin, M. Denhgan. Bernstein Ritz-Galerkin method 
for solving an initial-boundary value problem that combines neumann and 
integral condition for the wave equation. Numerical Methods for Partial 
Differential Equations, 26:1236-1246, 2009. 



20 



[4] E. H. Doha, A. H. Bhrawy, M. A. Sakcr. On the derivatives of Bernstein 
polynomials: An application for the solution of high-even-order deifferential 
equations. Boundary Value Problems, 2011. 

[5] J. Boyd. Chebyshev and Fourier Spectral Methods. Dover, 2nd edition, 2001. 

[6] B. Fornbcrg. A Practical Guide to Pseudospectral Methods. Cambridge 
University Press, 1998. 

[7] M. Y.Hussaini, C. L.Strcctt, T. A.Zang. Spectral Methods for Partial Dif- 
ferential Equations. NASA Contractor Report, NASA-CR-172248, 1983. 

[8] R. T. Farouki. Legendre-Bernstein basis transformations. J. Comput. Appl. 
Math. 119: 145-160, 2000. 

[9] http://code.google.eom/p/bernstein-poly/ 

[10] Y. C.Hon, W. Chen. Boundary knot method for 2D and 3D Hclmholtz and 
convection-diffusion problems under complicated geometry. International 
Journal of Numerical Methods in Engineering, 56:1931-1948, 2003. 

[11] M. Arad, A. Yakhot, G. Ben-Dor. A Highly Accurate Numerical Solu- 
tion of a Biharmonic Equation. Numerical Methods for Partial Differential 
Equations, 13:375-391, 1997. 

[12] D. A.Kopriva, J. H.Kolias. A conservative staggered-grid Chebyshev mul- 
tidomain method for compressible flows. Journal of Computational Physics, 
125(1):244-261, 1996. 



21 



